Self-sustaining oscillations of a falling sphere through Johnson-Segalman fluids 
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We confirm numerically that the Johnson-Segalman model is able to reproduce the continual 
oscillations of the falling sphere observed in some viscoelastic models, [HE]. The empirical choice 
of parameters used in the Johnson-Segalman model is from the ones that show the non-monotone 
stress-strain relation of the steady shear flows of the model. The carefully chosen parameters yield 
continual, self-sustaining, (ir)regular and periodic oscillations of the speed for the falling sphere 
through the Johnson-Segalman fluids. Furthermore, our simulations reproduce the phenomena 
characterized by Mollinger et al. pQ : the falling sphere settles slower and slower until a certain point 
at which the sphere suddenly accelerates and this pattern is repeated continually. 



I. INTRODUCTION 



Complex fluids can afford a variety of new phenomena, among others, the flow instabilities due to the internal 
structures that induce nontrivial interactions between fluids and macromolecules therein 3 , 4 . Typical flow instabil- 
ities observed in viscoelastic flows can be categorized into two distinguishing types: elastic instability and material 
instability 3J. Whereas the elastic instability in the polymer fluids is due to the strong nonlinear mechanical proper- 
ties of the polymer solutions [5] , the material instability is believed to have its origin in non-monotone dependence of 
stress on strain rate [3] . Examples of material instability, such as the shear banding [6j [7] , are observed in wormlike 
micellar fluids, that are fluids made up of elongated and semiflexible aggregates resulting from the self-assembly of 
surfactant molecules |S]. 

A recent and striking evidence of the flow instability in the wormlikc micellar fluids is discovered by Jayaraman 
et al., 2J: a sphere falling in a wormlike micellar fluid undergoes continual oscillations without reaching a terminal 
velocity; see also [9 . The continual oscillations of a falling sphere in a wormlike micellar fluid is quite exotic and in 
contrast to the case of generic polymeric fluids jTUHH] . Note that the oscillation of a falling sphere has been observed 
as well in another systems by Mollinger et al., [1] and was stated as an unexpected phenomenon. The steady shear 
flow of the wormlike micellar fluids, in which the falling sphere continually oscillates, is observed to display a flat 
region FIG. 2(a)] in stress-strain relation, which is known to be linked with the non-monotone shear stress-strain 
rate relations. 

The non-monotone stress-strain relations have been extensively investigated for such instabilities as shear banding 
[TUTS], shark-skin and spurt [THl E] and the Johnson-Segalman (JS) model [TS] showing such nonmonotone stress 
dependence on the strain has successfully provided qualitative agreements with experimental results in these instances 
[15] . It is then conjectured in [2] that "the continual oscillations of the falling sphere could be due to the same 
instability relevant to the non-monotone stress-strain relations" and suggested the JS model may produce the falling 
sphere experimental results in the wormlike micellar fluids. It is, however, apparent that the velocity fields of a falling 
sphere in wormlike micelles are much more complicated than the steady shear flows; the fluids are both elongated and 
sheared in the wake and near the sphere surface. 

In this paper, through a sophisticated numerical technique, we demonstrate that the minimal ingredient that 
contributes to continual oscillations is related to the non-monotone stress dependence on the strain rate. Moreover, 
we provide numerical evidences that the flow instability observed in falling sphere experiment [5] can be attributed 
to the material instability manifested by such a plateau region displayed in the stress-strain rate curve. The rest of 
the paper is organized as follows. The model description and a stable numerical scheme are presented in Section |H| 
A number of numerical experiments as well as their physical relevance are presented in Section |III| Finally, some 



concluding remarks are given in Section IV 



II. FLOW MODELS AND NUMERICAL SCHEMES 



We assume that a solid sphere of radius r s with density p s is accelerating with speed U s from rest under the influence 
of gravity g in a viscoelastic fluid of density pf with zero shear viscosity n contained in a cylinder of radius r c . Let 
p = Ps/pf be the density ratio, the symbols D t and d t denote the material time derivative and the time derivative, 
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FIG. 1: (Left) The plot of the shear stress as a function of the strain rate for different slippage parameter § = 1 — a. (Right) The plot 
of the shear rate at which the local maximum of the shear stress occurs (square) and the local maximum of the shear stress (circle) as a 
function of the slippage parameter £. Parameters used are Re = 0, fi s = 0.03, Wi = 0.45 and £ = 0.1 ~ 0.8. 



respectively. And, e z is the unit axial vector. The dimensionless momentum equations can be given as follows: 

-Vp + V ■ t + Re(d t U) e z 



ReD t u 
V ■ u 



0, 



(1) 



where U is the dimensionless speed of falling sphere, the parameters Re = PfV^r s jT] and p, s — r] s /r) and the total 
stress t — /i s (Vu + Vu T ) + t p . For characteristic scales, we used the following characteristic scales [13]: (1) length 
scale: radius of the sphere r s , (2) velocity scale: Stokes terminal velocity LTjv = 2r^(p s — pf)g/9r/K with r) = ij s + rj p , 
(3) time scale: r s /UN, and (4) stress and pressure scale: riUN/r s . 

The dimensionless JS constitutive equation for the elastic stress r p reads 



(2) 



T. p + Wi S e t p = m p (Vu + Vu T ), 
where Wi = AU~N/r s , M P = Vp/Vi an d the Gordon-Schowalter convected time derivative [3] is defined as 

5 e t p := D t T p - U(u)t p - T p ll(u.) T (3) 
with lZ(u) = |((a + l)Vu+ (a — l)Vu T ). The parameter a is related to the slippage parameter ( = l- aG [0, 1] [!§] ■ 



(b) 





FIG. 2: A schematic of the domain of the cylinder of radius r c in which the solid sphere of radius r s falls with the speed U(t) at time t. 

Note that, when £ = 0, the JS model reduces to the Oldroyd-B (OB) model [? ]. 
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FIG. 3: Speed of the falling sphere through the Johnson-Segalman models in two different meshes (the fine mesh is refined from the coarse 
mesh by the regular refinement). (Left) Parameters used are Re = 0.0325, /is = 0.59, Wi = 0.5, a = 6.3 and a = 0.8 and the difference 
between terminal velocities is 0.002538 (Right) (a) The coarse mesh, Mesh I for the computational domain with the aspect ratio p = 4.115 
in Figure [2] The fine mesh, Mesh II obtained from Mesh I. (c) Close up view of the mesh for domain with the aspect ratio p = 2 near the 
sphere, (d) Close up view of the mesh for domain with the aspect ratio p = 2 of Mesh II near the sphere. 

The equation of motion ([I]) is formulated in the moving frame viewed from the sphere [13] and is coupled with the 
dynamic force balance equation of the falling sphere, whose dimensionless form is given by (2Rep/3)dtU — 3K + Fd 
where Fd is the drag force [20] and K is the wall correction factor. The constitutive equation ^ can be reformulated 
in terms of the conformation tensor c [2T] , the ensemble average of the dyadic product of the end-to-end vector for the 
dumbbells as c + \N\5ec = (/i p /aWi)5, where S is the identity tensor. The conformation tensor c of JS model can be 
shown to be positive definite at all time |22H24j . In our numerical simulations, we assume that the flow is axisymmetric 
and apply the numerical methods proposed in [22], which is based on the semi-Lagrangian method (SLM) [25j [26] 
and preserves the positivity of the conformation tensor in the discrete level. Let (u old ,p old , c old ) denote the previous 
semi-discrete solutions. The backward Euler method combined with the SLM method applied for both D t u and D t c 
result in the following semi-discrete systems which will define the current time level solutions (u now ,p ncw , c new ). 

(Re/h t )u ncw + Vp ncw - ^ s Au ncw = (Re/h t )u old o y + V • c old + Re(d t U) old e z , (4) 

and Ac DCW + c ncvj A T = C, with A = 1/2 + Wi V 2 - ^(u llcw ) and C = (fi, p /a\N\)S + (W\/h t )c olA o y, where y denote 
the characteristic feet of the particle X at the current time, which can be obtained by solving the flow map equation 
that d s y(X,s) = u(y(X,s),s) and y(X,t) = X, fit is the time step size, (dtU)° ld is the acceleration of the falling 
sphere at i old . 

Note that the semi-discrete constitutive equation for the constitutive equation can be viewed as the Lyapunov 
equation |27) . For the spatial discretization, we apply the standard finite element approximations: continuous piecewise 
quadratic polynomial approximation for the velocity field, continuous piecewise linear approximation for the pressure, 
and continuous piecewise linear polynomial approximation for the conformation tensor. More detailed description of 
our algorithms can also be found at the review paper [55] ■ Note that the choice of our finite elements for the velocity 
and pressure has been shown to be stable for the axisymmetric Stokes equations [29| recently. For all simulations, we 
use the time step size h t = 10 -3 . The numerical solutions have been shown to converge with respect to the mesh size 
in case they achieve the steady-state; FIG. [3] (Left) indicates this visibly and the mesh convergence is observed. 

III. NUMERICAL EXPERIMENTS AND DISCUSSIONS 

For the numerical experiments, we focus on the choice of slippage parameter denoted by £. In particular, we have 
chosen parameter ranges that show the non-monotone stress-strain rate for the one dimensional parallel shear flows, 
for which the shear stress and strain rate k relation is t(k) = /j, s k + (£/jk)/[9£ 2 + £(2 — £,)k 2 }, where C = 1/Wi and 
(j, = (1 — /x s )/Wi. We note that the non-monotone stress-strain relation can always be obtained for £ > 0; see Fig. [l] 
(Left). It is observed that not all such parameters yield continual oscillations. For example, while £ = 0.2 gives the 
non-monotone shear stress-strain relation, it does not yield continual oscillations, which indicates that the choice of 
parameters based on the simple shear flows can be misleading. Our results, however, demonstrate these curves are 
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FIG. 4: The contour plot of velocity fields u (Left), v (Middle) and pressure, p (Right) from the computation of the Johnson-Segalman 
model at an instant time level, t = 25.665. The parameters are Re = 0.0325, /j. s = 0.03, Wi = 0.45 and the density ratio a = 6.3. 




FIG. 5: The contour plot of gradient of velocity fields d r u, d z u, d r v and d z v viewed near from sphere from left to right, respectively of 
the Johnson-Segalman fluid at time level t = 8.5367. Here (u,v) are the radial and axial velocity component of u. The parameters are Re 
= 0.0325, iJ, s = 0.03, Wi = 0.45, p = 6.3 and a = 4.115 

apparently related to the observed flow instability in the wormlike micellar fluids. More precisely, by varying the 
parameter £ = 0.0 ~ 0.8, we observe that the larger £ is, the smaller the local maximum value of the shear stress at 
which the transition to the negative slope of the shear stress curve occurs; see FIG. [T] 

The range £ = 0.4 ~ 0.8 produces continual oscillations and the amplitude of the oscillations is a increasing function 
w.r.t. £; see FIG. [8j The parallel shear flows have been studied by Renardy [T7], which indicate that the apparent 
local shear banding observed in FIG. [5] can in fact, be related to the flow instability. Unlike one dimensional shear 
flows [30j [31] , it has been shown that the shear flows is of the short-wave instability and unstable with respect to the 
two dimensional perturbation. Our numerical solutions further support important aspects of theoretical results on 
the material instability by Renardy [T7] , namely, the locality of the material instability. The instability in the parallel 
shear flow of JS model is shown to possess the feature: "the most unstable mode is found at the interfacial mode 
and the short-wave instability decays exponentially fast away from the interface position and the effect is localized 
in a boundary layer near the interface, not disturbing the bulk of the flow" . We consider two points: Xi is located 
near the sphere, whose distance from the surface of the sphere is 0.2293 and x^ is chosen away from the sphere, near 
the top wall and compare the solution behavior by plotting sample variables, the velocity fields u = (it, v) and the 
pressure p in the FIG. [6| Whereas at x\, all solutions are unsteady, at X2 the velocity component u is steady. This 
fact may justify our assumption on the flow symmetry. 

Our numerical results agree qualitatively with the experimental observations, [TJ [5] that the large sphere oscillates 
whereas the small sphere does not oscillate, which indicates that the continual oscillations may be affected by the 
wall effects. More precisely, we observed that the class of parameters, which produced the continual oscillation of 
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FIG. 6: (Left) The velocity fields u (solid), v (dot) and pressure p (solid dot) as a function of time at specified point x\ in the domain. 
(Right) The plot of u— component of the velocity u at the point X2- The parameters used are : Re = 0.0325, /j, s = 0.03, Wi = 0.45, a = 0.3 
and the aspect ratio is a = 4.115. 




falling sphere for the aspect ratio a = 4.115, did not produce the continual oscillations for the aspect ratio a = 8.115 
while the amplitude of the oscillation for a = 6.115 is smaller than that for the case a = 4.115. Furthermore, the 
oscillation pattern is clearly similar to the physical observations stated by Mollinger et al., pQ. Namely, the falling 
sphere settles slower and slower until a certain point at which the particle suddenly accelerates and this pattern is 
repeated continually, see FIG [7] and FIG [3] Note that this characteristic can be found at the work by Jayaraman et 
al., [2] as well. 

The flow past a sphere in the wormlike micellar fluids is observed to generate negative wake repeatedly before the 
falling sphere accelerates during the oscillations [2], i.e., the fluid behind the sphere moves in the opposite direction 
to the falling sphere, [101132]. The JS model produced the temporarily developed negative wake which disappears in 
time, see FIG. [9j This observation is similar to what is observed for the time dependent simulations of OB model, 
[55] and it may be not surprising since the formation of the negative wakes are observed for the class of FENE 
models [55J [51] at sufficiently low values of the extensibility parameter L of the dumbbell [55J [51] and both OB and 
JS models do not impose the finite extensibility of the dumbbells. This result indicates that the negative wake does 
not appear to be closely related to the continual oscillations of the falling sphere. Finally, unlike what is observed in 
physical experiments, that oscillations can occur for larger density ratios between the sphere and the fluids; while 
a Delrin sphere of density p s — 1.35<?/cm 3 with the diameter d = 3/16 in does not oscillate, a Teflon sphere of density 
p s = 2.17g/cm 3 of the same diameter does oscillate. As demonstrated in FIG.[9j the results shows that the simulation 
of JS model with three different density ratios, a = 1.3,4.3 and 8.3 produce similar pattern, which provide a room 
yet to be improved in mathematical modeling of the wormlike micellar fluids. 
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FIG. 9: (Left) Normalized speed of the falling ball in JS fluids as a function of time varying the density ratio p from 1.3 to 8.3. (Right) 
The plot of fluid velocity along the cylinder axis for JS model in time evolution with density ratio p = 6.3. For both computations, the 
parameters used are : Re = 0.0325, /i s = 0.03, Wi = 0.45, a = 0.3 and the aspect ratio is a = 4.115. 



IV. CONCLUSIONS 



In conclusions, the JS model is shown to reproduce the continual oscillations at certain range of parameters showing 
the non-monotone stress-strain rate relations, which is, therefore, believed to be the minimal ingredient to trigger the 
flow instability relevant to continual oscillations. Some qualitative pattern of the continual oscillation of the falling 
sphere has been observed as well. For more quantitative agreements with the real experiments, one has to employ the 
fundamental modeling and simulations of the microstructures of the wormlike micelles such as breaking and reforming 
of the wormlike micelles. This will be the worthy challenges in the future research. At this point, we are performing 
our simulation in the full three dimensional setting in order to eliminate the assumption that the flow is symmetric. 
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